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Abstract 

The aim of this paper is to contribute to the basic understanding of neuronal po- 
larization mechanisms by developing and studying a reaction-diffusion model for protein 
activation and inactivation. In particular we focus on a feedback loop between PIS kinase 
and specific GTPases, and study its behaviour in dependence of neurite lengths. 

We find that if an ultrasensitive activation is included, the model can produce polar- 
ization at a critical length as observed in experiments. Symmetry breaking to polarization 
in the longer neurite is found only if active transport of a factor, in our case active PIS 
kinase, is included into the model. 

1 Introduction 

Neuronal polarization is a key process in brain development, whose molecular origin is still 
far from being completely understood. Clearly, the differentiation of neurites into axons (in 
general only one) and several dendrites is a key mechanism for the structure and function 
of the brain, and thus the process is likely steered by external cues. On the other hand 
various lab experiments with cultured hippocampal neurons (cf. [5l [9]) show that a robust 
polarization and differentiation also occurs without external cues, driven by the activity of 
GTPases and other proteins. We refer to ^ El O |25l |28] for detailed discussions of neuronal 
polarity. 

During early development a neuron first extends several unspecified neurites that show 
random episodes of growth and retraction resulting in a constant average length. At a later 
stage a polarization in the distribution of several proteins occurs and usually one of the 
neurites is specified as the axon, which results in the stabilization of microtubuli (cf. [34j) 
and a phase of fast growth, while the other neurites grow slower and later become dendrites. 
The molecular origin of neuronal polarity is a spatial reorganization of proteins towards the tip 
of the longest neurite, which seems to appear when one of the neurites has reached a critical 
length (depending however on the overall configuration). In particular the sequential activity 
of the GTPases RaplB, Cdc42, and Rho/Rac GTPases has been found to be of fundamental 
importance for the differentiation of an axon. 
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The aim of this paper is to gain further understanding of neuronal polarization and 
symmetry-breaking as observed in cultured neurons by mathematical modelling of the GT- 
Pase signalling cascades, diffusion, and possibly transport. In particular we want to tackle 
the following questions: 

• Can a feedback loop between GTPases and PI3 kinase that has been found in experi- 
ments lead to an inherent length-dependent polarization mechanism ? 

• Is a suggested ultrasensive reaction essential for the polarization ? 

• What can cause symmetry breaking and determine the strongly preferred polarization 
of signaling molecules in the longest neurite ? 

We want to understand these issues by modelling activation of PIS kinase and GTPases 
such as Rho, Rac, Ras, Cdc42, RaplB, which are key players in the establishment of polarity 
and the growth of axons (cf. [TOl [71 |26l |3T1 |35j). Similar pathways of GTPases are found 
also in the polarization of other cell types (cf. [23l |6]). In particular the estabhshment of 
polarity in eukaryotic cells such as yeast, which appears to be a reasonably simple model 
system, has recently attracted much attention in the applied mathematics community (cf. 
[22ll2Zl[l2llHl[l3l[14]), and significant progress has been made very recently in understanding 
basic mechanisms as well as in the mathematical analysis of such models (cf. [HI [181 [201 [19]). 
Modelhng the polarization of axons is a topic less studied (cf. [24]), probably also due to 
the higher complexity in neurons. On the other hand the length-dependent mechanism and 
the preferential polarization of signaling cascades to the longer neurites is a robust process 
that a model should be able to explain at least qualitatively. Moreover, the changing length 
provides additional variability that allows to rule out certain models. 

The key assumption of our modeling approach, which seems to be in good agreement with 
experiments (cf. e.g. [9J) is that none of the neurites is predetermined to become the axon, 
and there is no polarization in neurite lengths without polarization in signaling proteins, 
more precisely GTPases. As described above, we consider the growth (and retraction) in the 
initial stage as random, and thus there can be some polarization in neurite length. However, 
without polarization in molecular distributions it seems quite unreasonable that a neurite at 
a certain length can become stable and subsequently differentiates to the axon. Thus, we ask 
whether mathematical modelling of GTPase activation can explain the stabilization of the 
longest neurite when it has reached a critical length by stochastic growth. 

As we want to focus on the polarization of GTPase distributions it is natural to set up a 
reaction-diffusion model to study spatial variation. Noticing that overall growth is slow, it is 
natural to assume that the reaction-diffusion process reaches equilibrium at fixed length, in 
particular if the equilibria are close also after an incremental growth step. The latter applies 
in particular for stationary concentrations being spatially homogeneous and thus we can argue 
that polarization will only appear if at a certain critical length the homogeneous stationary 
solution becomes unstable and a polarized solution showing a peak on one end appears. 
There are several suggestions of extracellular and stochastic processes that influence polarity 
(cf. [21 131 [28]), which can be understood as small perturbations of the homogeneous state. 
Thus, we consider the modelling of GTPase activation on fixed domains, but investigate the 
influence of different domain sizes. In particular we expect the homogeneous stationary state 
to become unstable at a critical length, which we interpret as length-dependent polarization. 
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The remainder of the paper is organized as fohows: In Section 2 we wih develop a basic 
mathematical model of activation mechanisms in developing neurites, which is further sim- 
plified in Section 3 and investigated numerically. In Section 4 we discuss the extension to an 
appropriate setup for a neuron with two neurites and investigate the possibility of symmetry- 
breaking. Finally we provide further analytical insight via a minimal model in Section 5 and 
conclude in Section 6. 

2 A Model for Activation of GTPases in Neurites 

In the following we derive a reaction-diffusion model for the polarization in neurites starting 
with a general approach to the modelling of GTPase activation and then proceeding to a 
specific signal pathway, which results in a system of reaction-diffusion equations. 

2.1 General Aspects in Modelling GTPases 

As mentioned above GTPases, appearing in cytosolic and membrane-bound forms, are of cen- 
tral importance for regulating neuronal polarization via feedback-loops and signaling growth. 
We therefore start by discussing the general modelling of GTPases first before turning to 
the specific signalling cascade. We follow the basic approach of [T?^ T3] and first detail a 
submodel for a single GTPase (ignoring for the moment the signal pathways between dif- 
ferent molecules), composed of three states the GTPase can switch to. There is an inactive 
cytosolic form (bound to GDI), an inactive membrane-bound form (bound to GDP), and an 
active membrane-bound form (bound to GTP). For the concentrations we use the following 
notations: 

• Ga denotes the membrane-bound active form 

• Gm denotes the membrane-bound inactive Form 

• Gc denotes the cytosolic inactive form 

We base the modelling on similar assumptions as |12| . motivated by experimental obser- 
vations: 

1. Neuronal polarization is sufficiently fast, also compared to growth. Thus we consider 
models on fixed geometries (however focusing on the impact on different neurite lengths) 
and assume that Rho-proteins are neither synthesized nor degradated at the relevant 
time scale. We only model the change between different forms. 

2. Each Rho-protein has spectific rates of activation and inactivation. These processes are 
indirect and controlled by GEFs and GAPs, respectively (see Figure [T]). 

3. The amount of GDI in the developing neuron is sufficiently high such that it does not 
become a limiting factor, and the change between the cytosolic and membrane-bound 
form is sufficiently fast. The latter is not crucial, but will allow to reduce the model to 
two forms, making the model easier to treat. 

4. Cytosolic forms diffuse much faster than membrane-bound forms (by a factor of around 
100, cf. ^3). Again this is not crucial for the basic modelling, but will become crucial 
for forming a Turing-type instability: the concentration of cytosolic forms will be rather 
homogeneous, while active forms will polarize at the membrane. 
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5. Since Rho-Proteins have almost equal sizes, we will assume that the diffusion coeffi- 
cients of different proteins in the active respectively inactive form are the same. This 
assumption is not crucial for the modelling and easily changed in the equations, but it 
appears convenient for numerical simulations. 

Let us denote the interior of the neurite as and the membrane as F C d^. The cytosolic 
concentration is modeled by a bulk diffusion equation 



in fi, 



(2.1) 



with A denoting the Laplace operator. The membrane-bound concentrations are modeled by 
reaction-diffusion equations on the surface 



dtGa — Dm^rGa + Ram{Ga, Gm) 

dtGm — Dm^vGm ~ Ram{Ga, Gm) + Rmc{Gm, Gc 



on r 
on r, 



(2.2) 
(2.3) 



where Ap denotes the Laplace-Beltrami operator on F. The reactions (activation and inacti- 
vation) between the two membrane-bound forms are given by 



RamiGa, Gm) = -k^GAPcGa + k+GEFoGm, 



(2.4) 



with the activation rate kQ and inactivation rate A:^, and the reactions (exchange) between 
the inactive forms are given by 

with the GDFcontrolled rates kQj^j und kQj^j (cf. Figure Finally, the boundary condition 
for Gc is given by 

DcVGc • n = -RmciGmi Gc) on F. (2.6) 



Membrane 




Cytosol 



Figure 1: GTPase reaction scheme (following ^6j). The GDFbound form diffuses in the cytosol 
and associates to the membrane with rate kQjj^. On the other hand the membrane-bound 
form is released and bound to GDI with rate 



GDI- 



We can further simplify the system based on the assumption of a fast exchange between 
the two inactive form, such that Kmc = on F, i.e. 



— ^GDI^^ ^GDI^<^ 



^ _ ^GDI ri 
'^GDI 



(2.7) 
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Finally, taking into account the elongated form of neurites, it is natural to reduce to a one- 
dimensional system. For this sake assume that x is the longitudinal axis of a neurite and let 
denote the two-dimensional plane where x = ^. Now we define effective concentrations of 
active and inactive forms 

Ua{x,t) = / Ga{x,y,z,t) ds{y,z), (2.8) 

Ui{x,t) = / Gm{x,y,z,t) ds{y,z) + / Gi{x,y,z,t)dydz. (2.9) 
JPxnr J JPxDQ 

Since we assume cross-sections to be small relative to the length, we can assume equilibration 
by diffusion in the y and z directions at a time scale relevant for longitudinal diffusion and 
reaction, such that we can well approximate all concentrations as independent of y and z. In 



particular we have with (2.7) 



ua - \Px n r|Ga, u, ^ (|p^ n n\ + n r\)Gm := ag^. (2.10) 

^GDI 

Thus, we find 

dtUa = / dtGa ds{y,z) ^ \Px^mDmdxxGa+Ram{Ga,Gm)) ^ DmdxxUa+Ram{Ua, X~^\Px^^Ui 



and 

dm 



PxHT 



dtGm ds{y, z) + / / dfGi dy dz 
J JPa-nn 



\Px n T\{DjndxxGm ~ Ram{Gai Gm)) + Dc\Px H fl\dxxGc 
Dfncdxx^i Rami'^ai A \Px ^ r|'U^), 



with the effective diffusion coefficient 

Dmc = x~^\Px nT\Dm + x-^\Px n n\Dc^^. (2.11) 

^GDI 

Thus, with redefining fcj = X~^\Px D r| A:^, we end up with the system 

dtUa = DmdxxUa- k^GAPcUa + k^GEFGUi, (2.12) 
dfUi = DmcdxxUi + k^^GAPcUa - k-^GFYcUi (2.13) 

for the activated and inactivated forms of the GTPases. Note that under typical parameter 
values we still expect D^c to be much larger than i?^ due to its dependence on Dc- 



2.2 Modelling the PI3-Kinase Signal Pathway 

In the following we try to model the most important signal pathways leading to neuronal 
polarization. We mainly follow signal cascades described in [21 [71 [351 EH]- Motivated by 
polarization models for eukaryotic cells, we particularly look for a feedback mechanism, which 
can lead to bistability or Turing-type instability. Such a feedback loop is known for the 
activation of PI3 kinases, which is discussed as a possible source for neuronal polarity [2, 7J. 
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During polarization PISkinases generate a feedback loop in a process involving their prod- 
uct PIP3, the GTPases Ras, RaplB, Cdc42, Rho, and Rac, as well as complexes such as 
PAR3/PAR6/aPKC. In order to keep the modelling at a reasonable complexity, we will not 
go into details of the biochemistry, but just assume some activity of canonical complexes 
involved in all the inhibitions and activations. The basic signalling pathway is sketched in 
Figure [2j where the Ki denote complexes involved and the indices a, i, c denote active, in- 
active membrane-bound, and cytosolic form, respectively. As demonstrated in the previous 
section, we will reduce all GTPases to an active and an inactive form. 

The GDP/GTP exchange of Rac is controlled by Cdc42 via the PAR3/PAR6/aPKC com- 
plex and recruitment of specific GEFs (like TIAM and STEF). We thus use a complex Ki 
for the activation of Rac in our model. PIS kinase is activated by the active forms of Ras 
and Rac, and inhibited by PTEN, which is controlled by Rho. We thus assume a complex 
K2, which contributes to inactivation of PIS kinase. Moreover, PISkinase indirectly activates 
Cdc42 via PIPS and RaplB (cf. [26j). In addition to those factors specific GAPs for Cdc42 
and Rac can intercept the feedback loop. We will denote them as complexes Ks und K4. 

We will model the above reactions simply via the law of mass action, except one, namely 
the activation of PIS kinase by Rac, for which we use a square law (alternatively Hill-type 
activation as e.g. in [19] are possible). Frequently ultrasensitive reactions, i.e. significantly 
steeper than the shape of Michaelis-Menten kinetics, are supposed. Indeed there is an experi- 
mental observation of ultrasensitivity of RaplB against signals of different receptors reported 
in [1^. Since the RaplB signalling pathways is related to PIS kinase, we assume its ultra- 
sensitive activation by Rac in our model. However, we think that also other ultrasensitive 
reactions in the feedback loop PIS kinase - Cdc42 - Rac would have a similar effect - a detailed 
determination of those remains a challenge for future experiments. On the other hand the 
existence of an ultrasensitive reaction seems crucial as we shall show in a minimal model in 
Section [5l 

The above considerations yield the following reaction-diffusion model for the PIS kinase 
concentration 



dtPISKa =k^RaSaPI3K^ + k^RaclPISK^ - k^K2PI3Ka + Dn^KAxPI^Ka, (2.14) 
dtPnK, = - k^RaSaPISK, - k^RaclPI3K, + k^K2PI3Ka + DpjskAxPI^K,. (2.15) 



Using the above reaction scheme in addition to the general submodel for GTPases from the 
previous section, we obtain a system of eight differential equations for the active and inactive 
forms of Ras, Cdc42, Rac, and Rho: 



dtRaSa ^k^^^Rasi - k^^^RaSa + DmdxxRaSa, 
dfRasi = - k^^^Rasi + k^^^RaSa + DmcdxxRasi 



(2.16) 
(2.17) 



dtCdca ^k^^^Cdci + k^PI3KaCdci - k^^^Cdca - k^K^Cdca 
- k^Cdca + k^Ki - k^Cdca + k^Ks^ DmdxxCdc, 

dtCdci = - k^^^Cdci - k^PI3KaCdci + k^^JJdca 
+ fcg K^Cdca + DjjicdxxCdci^ 



(2.18) 



(2.19) 
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Figure 2: Signalling scheme of GTPases involved in our model including the feedback loop 
PI3kinase - Cdc42 - Rac. 



dfRaCa ^k^^^Raci + k'l^KiRaci - k^ac^aca - k^^K^RaCa 

- k^RaCa + k^K4 + DmdxxRaCa, 
dtRaci = - k^^^Raci - k^^KiRaci + k^^^RaCa + k^^K^RaCa 

~1~ I-^mc^xxR^^i') 



(2.20) 
(2.21) 



dtRhoa ^k^ho^hoi - kj^^^Rhoa - k^Rhoa + k^ K2 + DmdxxRhoa, (2.22) 
dtRhoi = - k^^^Rhoi + k'^^^Rhoa + DmcdxxRhoi, (2.23) 

Finally we obtain model equations for the different complexes, for which we ignore diffusion: 

dfKi ^k^Cdca - k^Ki, dtK2 = k^Rhoa - fc2^i^2, 
dtK^ ^k^Cdca - k^Ks, dfK^ = k^RaCa - k^K^. 

We mention that the model can be seen as a mass-conserved reaction-diffusion system sim- 
ilar to those proposed for eukaryotic cells. In particular the following conservation properties 
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hold in the absence of diffusion: 



dt{RaSa + Rasi) = 0, (2.24) 

dtiCdca + Cdci + Ki+ Ks) = 0, (2.25) 

dt{RaCa + Raci + K4) = 0, (2.26) 

dt{Rhoa + Rhoi + K2) = 0, (2.27) 

dt{PI3Ka + PI3Ki) = 0. (2.28) 



We mention that the result of our modelling is a system of 14 equations, which is still difficult 
to handle. We will therefore try to further reduce it to a system of six equations for the 
active and inactive forms of PIS kinase, Cdc42, and Rac, whose feedback loop shall be the 
key aspect for polarization. 



3 A Six-Equations Model 

In the following we shall further reduce the system. First of all we assume that the reactions 
via all complexes are fast at the relevant time scales, which yields 

Ki =-^Cdca := hiiCdCa, K2 -^Rhoa := K.2RhOa, 

k^ 

Ks =-^CdCa := KsCdCa^ K4 = -^RaCa := n^RaCa. 

Moreover, we try to eliminate the equations for the GTPases Ras and Rho, since they act 
on PIS kinase, but are not influenced by the components of the feedback loop - we therefore 
assume they equilibrate and obtain the relations 

RaSa = -^Rasi (3.1) 

^Ras 

Rhoa = -^Rhoi. (3.2) 

^Rho 

Further approximating the inactive concentration as homogeneous and using the conservation 
properties we find 

Rasa = -^^^Ras^ := -fiRas^, (3.3) 



1 + 



^Ras 
'^Ras 



'^Ras 
^Rho 



Rhoa = T^Rho^ •= 72Rho^, (3.4) 

-1 I I^T^^c I ^E?^o 



1 + 7^ + 1^2 



^Ras 

^Rho ^Rho 



where the index * denotes the total concentrations in the system, which we assume fixed. 
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Inserting those simplifications into (2.14)-(2.24) we finally obtain a system of six reaction- 



diffusion equations for the concentrations of PIS kinase and the GTPases Cdc42 and Rac: 

dtCdca -k^dcCdc^ + k^PnKaCdc, - k-^,Cdca (3.5) 

- kg/^4RaCaCdCa + DmdxxCdca, 

dtCdc, = - k^dcCdc^ - k^PISKaCdc, + kcdcCdca (3.6) 

+ kg hi/^RaCaCdca + DjjicdxxCdci^ 

dtRaCa ^k^^^Raci + k^^i^iC dcaRaci - k^^^RaCa (3.7) 

- k^^nsCdcaRaCa + DmdxxRaCa, 

dfRaci = - k^^^Raci - k^^hiiC dcaRaci + k^^^RaCa (3.8) 

+ k^^KsCdcaRaCa + DmcdxxRaci, 
dfPISKa = - k^n27iRho^PI3Ka + k6-f2Ras^PI3Ki (3.9) 

+ krRaclPISKa + DpisKadxxPI^Ka, 
dtPISKi = + k^i^2liRho^PnKa - kGj2Ras^PI3Ki (3.10) 

- krRaclPISK, + DpjskAxPI^K,. 

3.1 Dimensionless Variables and Scaling 

In order to eliminate the dependence on dimensions and in particular to rescale the domain 
to unit size, which simplifies the systematic investigation of length-dependent polarization, 



we perform a scaling of the variables in the system (3.5)-(3.10). We denote the dimensional 
variables used in the previous section by x and t and introduce new space and time variables 
and new concentrations via 



X_ ^ _ Lq 

L' ^ ~ L 



ui — CdcaL^ Us — RaCaL^ — PISKaL 

U2 — CdciL^ u/^ — RaciL, uq = PI3KiL 

where r is a suitable time scale and L is the total length of the system (to be varied in the 
following). Lq ia an arbitrary but fixed reference length used to introduce the dimensionless 
parameter e and thus formulate the question of length- dependent polarization in terms of e. 
For e large (small length) we do not expect polarity, thus there should be a homogeneous 
solution u = {ui,U2,us,U4,u^,uq). At critical e, the homogeneous solutions should become 
unstable and different solutions showing maxima of the active concentrations {ui for i odd) 
at the boundary should appear. 



The scaled version of ( |3.5| )-( |3.10[ ) is given by 

dfUi —aiU2 + ea2U^U2 — a^ui — ea^u^ui + die^dxxUi, (3.11) 

dtU2 = - aiU2 - ea2U^U2 + a^ui + ea^usui + e^dxxU2, (3.12) 

dtUs =0:5^/4 + eaQUiU4 — a^us — ea^uius + die^dxxUs^ (3.13) 

dtU4 = — 0^5^x4 — eaQUiU4 + 0^7^x3 + ea^uius + e^dxx^^, (3-14) 

dtU5 = - cxqu^ + aioUQ + e^aiiuluQ + d2e^dxxU5i (3.15) 

dtUQ =0:9^5 - c^iqUg - e^aiiuluQ + dse^d^xUe, (3.16) 

with the new parameters given in Table [ij 
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ai = 
0^5 = k 



Cdc-^o/^ri 



Rac 

T = Ll/Dm 



«2 k^Lo/Dmc 
Q/4 k^K^Lo/Dm 

ae = k^QhiiLo/Dmc 
as = k^^i^sLo/Dmc 
«io = kQ-f2Ras^Ll/Dr^ 

ds DpjsKi I Dmc 



Table 1: Overview of dimensionsless parameters in (3.11 )-(3.16) 



We use this model now on the unit interval with no-flux boundary conditions, i.e. 

a,«iUo,i = 0. (3.17) 
Note that the system has three volume-conserving subsystems, we have 

dt / {ui -\- Ui^i)dx = 0, fori = 1,3, 5. (3.18) 
Jo 

Of course it would be an ideal case to obtain the parameters in the model from experi- 
mental observations, which are currently not available in most cases however. Based on data 
on protein diffusion used in [22j we use a factor 100 in the diffusion coefficients between the 
cytosolic and membrane-bound forms of GTPases. Their values for cytosolic diffusion coeffi- 
cients are around 10 — 50/xm^/5~^ and for membrane diffusion coefficients around O.lfim? /s~^, 
which we use for our simulations. For PI3 kinase we use the value of O.l/xm^/s"^, since it 
is still bound to the membrane via receptors. Moreover we assume cytosolic diffusion for 
the inactive form and thus a diffusion speed of lO/xm/s"^, motivated by its association with 
inactive RaplB. Rates of activation and inactivation cannot be found in literature, so we use 
rough estimates leading to reasonable time scales. The used dimensionless values are given 
in Table d 



3.2 Numerical Results 



In the following we present results of numerical simulations of the system (3.11 )-(3.16) on 



the unit interval with boundary conditions (3.17). All results were computed in MATLAB 



with the pdepe tool. To check that there are no discretization artefacts and that computed 
equilibria are robust we also used an independent upwind finite difference implementation 
with backward Euler time discretization, which led to the same results up to numerical errors 
and are not displayed here. 

The initial value in all simulations is a perturbation of homogeneous states for all con- 
centrations. Both cosine-type and random perturbations were used, which lead to the same 
results except that the direction of polarity is predefined by the sign of the cosine perturbation, 
while random perturbations lead to random choice of direction. 

In order to investigate the ability to reproduce length- dependent polarization we run the 
system with initial concentrations Ui = 0.5 + 0.001 cos ttx for different sizes of L until a 
stationary state is reached {t = 100000 seems to be a reasonable choice). The results are 
given in Figure [3) One observes that for the parameter set given in Table [2] length-dependent 
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Parameter Value Parameter Value 



ai 


0.001 


as 


0.01 


a2 


0.5 


ag 


0.1 


as 


0.01 


aio 


5 


Q!4 


0.01 


an 


0.001 


0:5 


0.001 


di 


0.01 


ae 


1 


^2 


0.01 


ar 


0.01 




0.01 



Table 2: Dimensionless parameter values used for numerical simulations 



polarization occurs at around L — 7. For smaller values the constant solution is stable, while 
for larger values a peak close to the left boundary (due to the higher perturbations on the 
left) appears, which is more and more pronounced for increasing L. To illustrate the time 
evolution leading to these stationary states we plot the concentrations at several intermediate 
time steps for the case L — 12. 



L=6.5,T=100000 



- Cdc42, 

- Cdc42i 



L=7,T=1 00000 



L=9, T=100000 L=12, T=100000 L=15, T=100000 

-, , , , , , 1 I , , , , , 1 o.14| 1 




0123456789 2 4 6 8 10 12 5 15 

X X X 



Figure 3: Stationary concentrations for different sizes of L. 



Clearly the specific perturbation costtx has a certain influence on the selection of the 
stationary solution, since e.g. for a perturbation with opposite sign one would end up with 
a solution mirrored at x = ^ due to the inherent symmetry in the system. However, apart 
from the selection of a solution with left or right peaks, the perturbation seems not to be too 
crucial. To test this we use various random perturbations (uniformly distributed in each grid 
point) and the resulting stationary solutions are always the same as for the low-frequency 
cosine perturbation, except that with same probability we also obtain the mirrored one. To 
illustrate this we plot the resulting stationary concentrations for L = 12 in the left panel of 
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Cdc42„ L=12 



Rac„ L=12 



PI3K„ L=12 




Figure 4: Evolution of active Cdc42, Rac and PIS kinase concentrations in a case leading to 
polarization. 



Figure [5| The only type of perturbation that does not yield peak solutions also for large L is 
a deterministic high-frequency one. This is not surprising since for large pure frequencies the 
diffusion dominates and smoothes the perturbation. This feature is illustrated in the right 
panel of Figure [5] with a plot of the resulting stationary concentrations in the case L — 12. 



L=12, T=100000 



- PI3K, 
PI3Ki 




L=12, T=1 00000 



- PI3K, 
RISK, 



Figure 5: Stationary concentrations for a (pointwise uniformly) random perturbation (left) 
and a high frequency perturbation (right). 



4 Symmetry Breaking and Axon Polarization 

In the previous section we have seen that the reaction diffusion model for PI3 kinase, Cdc42 
and Rac is indeed able to reproduce length-dependent polarization, which is a first confirma- 
tive answer to the possible existence of an inherent mechanism. The remaining question is 
how the symmetry in the system is broken, in particular how the observation that polariza- 
tion to a peak appears usually at the tip of the longest neurite (cf. [5]) can be explained. 
For this sake we extend the model in the following to deal with a setup including neurites 
of different lengths. Since we want to keep the setting reasonably simple we use a system 
with two neurites separated by the soma as sketched in Q, which still seems reasonable to 
be reduced to one spatial dimension and for which there is also some biological motivation as 
highlighted in [4J. 



12 




Figure 6: Setup of a neuronal cell with two developing neurons of lengths Li and L2, with a 
soma of length Lq inbetween. 



4.1 Modelling the Soma 

Our main rationale for including the soma into our model is that the mobility (and hence the 
diffusion coefficient) is lowered. For the cytosolic form this is due to the cell nucleus, which 
prevents flow in a large fraction of the soma, and for the membrane-bound form this is due to 
the additional curvature of the membrane along the soma, which is otherwise not incorporated 
in our model. Consequently we change all terms dxx^i in ( |3.11| )-(3.16) to dx{di{x)dxUi)^ with 




if < X < /i and (/i + 1^) < x <l, ^ ^ 

if /i<X<(/i+/o). 

For illustrative numerical tests we choose the values q = 0.6 for i odd (membrane) and 
Q = 0.3 (cytosol) and all other parameters as in the previous section, but similar results have 
been found for other settings. For more intuitive visualization we rescale the lengths on the 
X-axis in all plots and center the system such that the mid point of the soma is located at 
X = 0, in addition we use a light coloring of the soma region. We again use initial values 
Ui{x^ 0) = 0.5 + 0.001 cosTTX. 

Figure [7| demonstrates that still length-dependent polarization, in this case at around a 
length L — 7.3 of the overall system, appears. The increase in the critical length compared 
to Figure |3] is due to the fact that the overall diffusion is smaller due to the soma. In these 
examples polarization occurs in the (longer) left neurite. However, there is no symmetry 
breaking towards polarization in the longer neurite, but the effect is caused by the positive 
initial perturbation in the left part. In a variety of simulations with different perturbations 
and length distributions we observed similar behaviour as in the model without soma, in 
particular a nonsymmetric initial perturbation most strongly influences the choice of the 
polarization. This is illustrated in Figure [8] with results varying the lengths Li and L2, while 
keeping all other parameters - in particular the shape of the perturbation - as above. A 
conclusion from these simulations is that the existence of a soma is not sufficient to obtain 
a symmetry-breaking towards polarization in the longer neurite as observed in experiments. 
We will thus consider a further extension in the next section. 



4.2 Active Transport 

Besides the mechanisms included in the model so far, it has been suggested that active 
anterograde transport has an important role in neuronal polarization (cf. 1^51 [30]). We will 
thus extend the model to include such a transport term. Our detailed motivation is a report 
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of active transport of Shootinl in [17], which has been identified to be a protein related to 
the localization of PIS kinase (cf. [29]). Thus, we modify the (scaled) equation for the active 
PIS kinase concentration to 

dtu^ = - aQU^ + e'^aioujuQ + auUQ + e'^d^ {d^{x)dxU^) (4.2) 

with 

{-V if < X < /i, 
V if (/i + /o) <^</2, (4.S) 
else. 

Since we are not aware of another experimental proof of active transport we keep the equations 
for the other concentrations unchanged here, but remark that our computational tests show 
similar results if active transport is added to any other equation. We also mention that here 
we use equal constant velocities in both neurites for simplicity. In [SOJ it was even suspected 
that the transport velocity for Shootin 1 is growing with the neurite length, which also does 
not change our qualitative results. 

The surprising result of using (even small) active transport is that symmetry breaking 
appears as in experiments and all observed effects are reproduced. The polarization now 
occurs mainly in the longer neurite even for positive perturbations in the shorter one. This is 
illustrated by the time evolutions of the active concentrations in Figure [9} again with initial 
value Ui{x^ 0) = 0.5 + 0.001 costtx. The positive initial perturbation in the shorter left neurite 
reverted towards a positive perturbation in the longer neurite and subsequently a similar 
evolution towards the stationary state as in the symmetric case. 

Some stationary solutions in the case of active transport, all obtained with the same initial 



value, are shown in Figure 10 The polarization appears in any case in the longer neurite. 
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L=12, Li=7, 1^=4, T=1 00000 



L=12, L,=4, U=7, T=100000 



L=12, L,=10.5, U=0.5, T=100000 




Figure 8: Stationary concentrations for constant length L and varying neurite lengths Li and 

L2. 



Another interesting observation is that stronger transport always leads to polarization 
in the longer neurite as we illustrate by Figure 11, Here at a small transport v — 0.001 the 
polarization occurs in the left neurite, which is a bit smaller. However, doubling the transport 
speed also yields polarization in the longer one despite the very small length difference and 
the positive perturbation in the shorter one. 

Finally, an experimental observation that is reproduced with active transport as well is 
the fact that polarization depends on the length configuration of the neurites, in particular 
their length difference in the case of two neurites. This is illustrated by the results in Figure 
12 , where the total system length L — 7.4 is kept fixed. However, in the case of two neurites 
of equal length this does not suffice to produce polarization, whereas in the case of one very 
large and one very small neurite the polarization clearly occurs. 



5 Analysis of a Minimal Model 

In order to gain some analytical insight into the polarization behaviour we study a minimal 
model in the following by reducing to the essential qualitative features, namely an active and 
an inactive form together with a nonlinear feedback loop in the activation. 

Denoting the active concentration by ui and the inactive by we study a system of the 
following form 



dtui ^e^dxxUi + p{eui)u2 - Sui, 
dtU2 ^De^dxxU2 - p{eui)u2 + 5ui, 



(5.1) 
(5.2) 



for X G (0, 1) and t > 0. Clearly 5 > and p : 



is a potentially nonlinear function 
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Cdc42„ L=12, L,=5, U=6 Rac„, L=12, L,=5, U=6 PI3K„ L=12, L,=5, U=6 




Figure 9: Evolution of active Cdc42, Rac and PIS kinase concentrations in a case leading to 
polarization. 



modelling the activation feedback loop. We normalize the mass such that 



/ {ui + U2) dx = 1. (5.3) 



We analyze the existence and stability of stationary solutions of ( |5.l| ), ( |5.2| ). Homogeneous 
stationary solutions satisfy 

= (5.4) 

p{eul){l-u{) = 5u{. (5.5) 

Stability of homogeneous solutions is given if the Jacobian of the reaction term 

F{ui,U2) = {p{^ui)u2 - 5ui, -p{eui)u2 + Sui) 

has eigenvalues with nonpositive real parts. The eigenvalues are given by Ai = und A2 = 
— 5 — p + ep'u2. Hence, we obtain the stability condition 

5 + p{eul) > ep\eul)u^2' (5-6) 

To check for Turing instability it suffices to consider the linearization of the reaction diffusion 
system for perturbations cos(A:7r), with a natural number fc, since those are the eigenfunctions 
of the second derivative with Neumann boundary conditions. For p. = k^Tv^ the eigenvalues of 
the linearizations can be computed as in the case = and we obtain the following condition 
for the appearance of a Turing instability: 



Dpe^p'ul > Dpe^S + Dp^e"^ + pe^p. (5.7) 

p{eul) 



6u* 

We insert U2 = / and simplify for D ^ 1 to 



^^^^^^ > 1 + (5.8) 

For small L respectively large e the term is dominating, thus there is no instability. For 
sufficiently large L the instability condition is satisfied if 

p{w)w > p{w) for all w > 0. (5.9) 
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Figure 10: Stationary concentrations for constant length L and varying neurite lengths Li 
and L2. 



Now we observe that (5.9) cannot be satisfied for linear />, since this would always imply the 
opposite inequality for all w. Hence, some ultrasensitivity is needed to make the homogeneous 
solution unstable at some length L. 

We also see that with a quadratic atcivation the system can produce a length-dependent 
Turing instability. Consider 

p{w) = a + /3w'^, (5.10) 
then p\w) = 2^e^w^ hence instability can occur if 

^-1 = ^^-1 = =^^>0. (5.11) 
This is true in particular if a = 0. 



6 Conclusions 

We have studied the mathematical modelling of neuronal polarization focusing on an exper- 
imentally observed feedback loop between PIS kinase, Rac and Cdc42 and investigated the 
model behaviour by numerical solutions in a setup with two neurites separated by the soma. 
Our main conclusions are the following: 

• A feedback loop in the activation of GTPases and PI3 kinase can lead to an inherent 
polarization mechanism at critical length. However 

• Some kind of ultrasensitivity is needed to produce instabilities leading to polarization, 
the latter do not occur with simple linear reactions. 
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L=12, L,=5.45, L2=5.55, v=0.001 , T=100000 L=12, L,=5.45, ^=5.55, v=0.002, T=100000 




_, 1 , — <- < 'I 0' ' "I ' ' ' ' J 

-4 -2 2 4 6 -4 -2 2 4 6 



Figure 11: Stationary concentrations when varying the transport speed. 




Figure 12: Stationary concentrations for constant length L and varying neurite lengths Li 
and L2. 



• Symmetry breaking as observed in experiments can be obtained by active anterograde 
transport. This is however not specific to the transport of a particular species, and - 
since the mathematical form for domain growth is analogous to the active transport 
term - could even occur due to mechanical effects. Note the experiments of [15j, which 
lead to polarization caused by pulling on a neurite. 

Our results suggest a variety of future questions to be solved. Of course our model is far 
from being quantitative and the need for extensions as well as calibrations of model parameters 
are to be clarified from experiments. The role of active transport during development is 
another issue that calls for further experimental investigations. 

From a mathematical point it is a challenging issue to characterize the existence of non- 
trivial steady states in cases of Turing instability, which so far has been successful in only few 
cases in literature (cf. [32l|33]), not for the mass-conserving model suggested here. An even 
more challenging and fundamental question is to understand the effect of a transport term in 
the equation and the detailed reasons for the symmetry breaking. 
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